function T = CTETorque(  x, epsilons )
%CTETORQUE Summary of this function goes here
%   Detailed explanation goes here
alpha0 = x(1); r1 = x(2); r2 = x(3); r3 = x(4); r4 = x(5); k2 = x(6); k3 = x(7); k4 = x(8);

%compute C4LM parameters in undeflected state
delta0 = sqrt(r1^2+r2^2-2*r1*r2*cos(alpha0));
eta0 = acos((r1^2+delta0^2-r2^2)/(2*r1*delta0));
xi0 = acos((r3^2+delta0^2-r4^2)/(2*r3*delta0));
lamda0 = acos((r4^2+delta0^2-r3^2)/(2*r4*delta0));

beta0 = xi0-eta0;
gama0 = pi-lamda0-eta0;

%calculate torques generated by C4LM with different deflection levels
alphas = alpha0+epsilons;
deltas = sqrt(r1^2+r2^2-2*r1*r2*cos(alphas));
etas = acos((r1^2+deltas.^2-r2^2)/(2*r1*deltas));
xis = acos((r3^2+deltas.^2-r4^2)/(2*r3*deltas));
lamdas = acos((r4^2+deltas.^2-r3^2)/(2*r4*deltas));

betas = xis-etas;
gamas = pi-lamdas-etas;

%Calculate Generalised Coordinates
psi2s = epsilons-betas+beta0;
psi3s = gamas-gama0-betas+beta0;
psi4s = gamas-gama0;

%Calculate troques
T2 = k2.*psi2s;
T3 = k3.*psi3s;
T4 = k4.*psi4s;

h32s = (r2*sin(gamas-alphas))./(r3*sin(betas-gamas));
h42s = (r2*sin(betas-alphas))./(r4*sin(betas-gamas));

T = T2-(T2+T3).*h32s+(T3+T4).*h42s;
end

